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Abstract 

A  first  order  analytical  approximation  of  the  tesseral  harmonic  resonance  perturbations  of  the  Keplerian 
elements  is  presented,  and  the  mean  elements  (the  Keplerian  elements  with  the  long  period  portions 
averaged  out)  will  also  be  given  in  closed  form.  The  results  of  a  numerical  test,  which  compares  the 
analytical  solution  against  a  numerical  integration  of  the  Lagrange  equations  of  motion,  will  be 
summarized,  and  the  implementation  of  the  solution  in  the  analytical  orbit  determination  routine  ANODE 
(at  Lincoln  Laboratory)  will  be  outlined. 
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AN  ANALYTICAL  TREATMENT  OF  RESONANCE  EFFECTS 

ON  SATELLITE  ORBITS 


1.  INTRODUCTION 

Normally  the  longitude  dependent  tesseral  harmonics  in  the  expansion  of  the  Earth’s  potential  produce 
only  short  period  perturbations  of  satellite  orbits,  and  these  are  small  in  comparison  to  the  dominant 
latitude  dependent  zonal  harmonic  Jj.  Typically  analytical  orbit  determination  routines  (such  as  ANODE 
at  Lincoln  Laboratory)  carry  the  short  period  effects  to  o(J-,),  and  so  the  tesseral  harmonic  effects  are 
ignored.  When  the  mean  motion  of  the  satellite  is  nearly  commensurate  with  the  rotation  rate  of  the  earth, 
however,  the  trajectory  of  the  satellite  repeats  itself  relative  to  the  earth  and  the  perturbations  due  to 
certain  ‘critical’  tesseral  harmonics  build  up  at  each  passage  in  the  same  spot.  In  this  case,  there  can  be 
important  long  period  effects  which  should  not  be  ignored.  In  this  report,  we  will  present  a  first  order 
approximation  of  the  isolated  resonance  effects  by  integrating  the  Lagrange  equations  for  the  general 
perturbations  directly.  This  solution  is  first  order  in  the  sense  that  the  elements  on  the  right  hand  side  of 
the  Lagrange  differential  equations  which  are  not  expressed  in  terms  of  time  explicitly  are  held  fixed  and 
constant.  The  corresponding  mean  elements  will  also  be  given  in  closed  form,  and  the  results  of  a  test 
which  compares  the  solution  presented  here  with  numerical  integration  of  the  equations  of  motion  will  be 
summarized. 

The  problem  of  resonance  effects  on  satellite  orbits  was  treated  by  a  number  of  authors  throughout  the 
1960s.  The  reader  is  encouraged  to  consult  Allan  [1],  Garfmkel  [2],  andGedeon  13].  Garfinkel  solved 
the  so-called  ‘ideal  resonance  problem',  which  was  applicable  to  a  variety  of  resonant  situations  and 
demonstrated  how  interesting  and  complicated  the  notion  of  resonance  can  be.  In  many  ways,  almost  all 
of  the  current  literature  on  resonance  can  be  traced  back  to  Allan's  paper  [1].  Gedeon  expanded  upon 
Allan’s  analysis  and  presented  a  general  study  of  resonance  effects  on  satellite  orbits  by  integrating  the 
equations  of  motion  numerically.  It  is  Allan's  work,  however,  which  influences  the  present  study  more 
than  the  others.  Indeed,  almost  all  of  the  work  in  the  literature  on  resonance  effects  on  artificial  satellite 
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orbits  consider  only  the  long  period  change  in  the  ‘broken-legged’  longitude  (i.e.,  the  window  of  where 
the  satellite  is).  The  changes  in  the  other  Keplenan  elements  are  also  important  in  generating  an  accurate 
approximate  ephemeris,  and  these  changes  are  influenced  directly  by  the  acceleration  of  this  longitude. 
Allan’s  suggestion  [1,  p.  1332]  of  expressing  the  longitude  as  a  function  of  time  using  the  Jacobi  elliptic 
functions  has  spawned  the  idea  in  this  report  to  exploit  that  representation  of  the  longitude  to  yield  a 
complete  analytical  solution  of  the  long  period  resonance  effects  on  each  of  the  Keplerian  elements. 

The  solution  given  by  Garfmkel  [2]  applied  to  the  artificial  satellite  problem  should  yield  the  same 
results  as  the  solution  outlined  in  this  report.  Recall  that  Garfinkel  applied  the  von  Zeipel  method  to  solve 
the  equations  of  motion.  Our  solution  will  be  obtained  by  integrating  the  equations  of  motion  directly  so 
that  each  of  the  Keplerian  elements  will  be  given  as  explicit  functions  of  time.  The  mean  elements  and 
long  periodic  corrections  can  be  readily  obtained  by  manipulating  the  solution  directly.  Thus,  we  feel  our 
method  is  more  simple  to  understand  and  to  implement  in  existing  analytical  satellite  orbit  determination 
routines  where  there  is  assumed  to  be  no  interaction  between  the  various  harmonics  in  the  geopotential. 
Our  solution  will  work  particularly  well  with  the  well-known  classes  of  satellites  effected  by  resonance 
such  as  the  near  circular  synchronous  and  half-synchronous  satellites.  It  is  for  these  objects  that  this  study 
is  intended,  and  'he  work  presented  here  is  not  meant  to  be  an  over-all  treatment  of  the  universal 
resonance  problem. 


It  is  easiest  to  describe  the  appropriate  resonance  problem  by  considering  Kaula’s  expression  of  the 
geopotential  disturbing  function  in  terms  of  the  Keplerian  elements  [4,  p.  37]: 

oo  /  /  oo 

X  Vtmpq  ' 

1=2  /7i=0  p=0 

where 

0“A  W'>G/m('V,„{cos,(W  ffor  e"en  • 

a\a/  n  L sin(y ,  J  for  l-m  odd 

and  ' q 

Vlmpq  =  0  ~  2p)v>  +  0-2 p  +  q)M  +  -  0  -  \,m). 

In  the  second  equation,  Flmp(i)  is  the  inclination  function.  Gtpq(e)  is  the  eccentricity  function,  p  is 

Newton’s  gravitational  constant  times  the  mass  of  the  earth,  +  sf  is  the  (unnormalized) 

coefficient  of  the  spherical  harmonic  of  degree  /  and  order  m,  Xlm  =  itan  ~](SlmlClm)  is  the  corresponding 
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reference  longitude  along  the  equator,  a  is  the  semimajor  axis  of  the  satellite  orbit,  and  ae  is  the  (mean) 
radius  of  the  earth.  In  the  equation  for  the  angular  argument  V|/,  M  denotes  the  mean  anomaly,  co  is  the 
argument  of  perigee.  0  is  the  mean  Greenwich  Sidereal  time,  and  Q  is  the  longitude  of  the  right  ascending 
node. 

We  will  write  s  =  s0  +  As  to  denote  the  mean  motion  in  revolutions  per  day,  where  s0  is  integral  and 
|Ax|  <  1/2.  Now  exact  orbital  commensurability  with  the  earth's  rotation  can  be  expressed  as 

M  +  to  -  jo(0  -  ft)  .  (1) 

Of  course,  we  are  using  the  dot  over  the  eiement  to  denote  differentiation  with  respect  to  time.  With  the 
above  equation  as  motivation,  we  introduce  the  following  fundamental  quantity  of  longitude  through  the 
differential  equation 

•  Ar  •  • 

X  =  —(0  -  Q)  . 
s0 

and  by  assuming  ( 1 )  we  write 

X  =  —  (M  +  co)  -  (0  -  Q.)  . 
s0 

Thus,  we  can  rewrite  tty  as 

VlmpqM  =  (^1  -  +  q  -  +  03)  +  m(X  -  Xlm)  -  qi 0  .  (2) 

The  underlying  critical  indices  will  be  those  sets  of  /,  m,  p,  q  which  satisfy 

I  - 2p +  q~  —  -  0  . 
s0 

The  quantity  X(t)  can  be  physically  interpreted  as  the  osculating  value  of  the  longitude  of  the  ascending 
node  of  the  mean  satellite  [3]. 

If  we  consider  only  the  critical  tesseral  harmonics  in  the  disturbing  function,  then  the  disturbing 
function  will  take  on  the  following  form: 

(3) 

crit. 

where 

Vlmpc  =  -  q«>  ■ 
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To  simplify  the  following  discussion,  we  will  temporarily  assume  *hat  tj limpq  =  mX\  that  is,  we  will 

assume  that  either  q  -  0  or  else  that  d)  =  0.  One  should  see  [3]  for  a  more  general  treatment.  At  any  rate, 

these  assumptions  can  safely  be  made  if  the  eccentricity  is  very  small  or  if  the  inclination  is  close  to  the 

critical  inclination  (cos/  =  ±  l/s/5).  Now  from  the  Lagrange  equation  of  the  general  perturbation  of  the 
•  2  3/? 

semimajor  axis  a  -  — — .  it  follows  from  (3)  that 
1  nadM 


,F>mPWlpq{e\Jlm{ 


-sin (\iimpq)  for  /  -  m  even 
COS(V/mM)  for  1  ~  m  odd 


This  change  in  the  semimajor  axis  will  in  turn  cause  an  acceleration  in  the  longitude  X 


X 


—  l~)i 

da\s0J 


,F,mPW<Pq^iJ-^lmPq)  for  l  -  m  even 
L  cos (\fimpq)  for  /  -  m  odd 


(4) 


We  have  assumed  that  the  accelerations  of  (0,  C2.  and  8  are  all  negligible.  In  an  attempt  to  analytically 
integrate  (4),  we  will  assume  that  the  elements  a,  e.  and  /  on  the  right  hand  side  remain  fixed  and  constant. 


Once  this  assumption  is  made,  a  first  integral  follows  easily 
crit.tf  -Tq 


(5) 


where  C  is  a  constant  of  integration. 


The  above  differential  equation  (5)  can  only  be  integrated  further  analytically  if  we  isolate  one  critical 
tesseral  harmonic  in  the  disturbing  function  (3).  If  this  assumption  is  made,  a  global  approximation  is 
subsequently  possible  under  the  assumption  that  no  two  critical  tesseral  harmonics  interact.  For  the 
remainder  of  this  report,  we  will  restrict  ourselves  to  this  situation. 
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2.  THE  SOLUTION  OF  THE  SIMPLE  PENDULUM  PROBLEM 


The  simple  pendulum  problem  in  physics  can  be  represented  mathematically  by  the  relation 

y  =  -£?2siny  ,  (6) 

where  Q  is  positive  and  constant.  The  solution  is  not  difficult  to  obtain,  and  we  will  outline  the  important 
details  for  reference. 

A  first  integral  is  found  as  in  (5)  and  is 
y2  =  C  +  202cosy  . 

The  integration  constant  C  can  be  expressed  in  terms  of  the  initial  conditions  y0  and  y0,  and  we  can  write 
C  =  y02  -  2£}2cosy0.  Clearly  C  >  -  2£52cosy,  and  so  if  C  is  small  enough,  then  y  will  be  prevented 
from  making  a  full  cycle  from  0  to  2ji.  If  we  write  cosy  =  1  -  2sin2  (^) .  then  the  sine  of  the  maximum 

deviation  from  the  focal  point  of  libration  (if  such  a  value  exists)  corresponds  to  y  =  0.  Explicitly,  we 
have 


.  i  V™ 

sm-  <— ) 


C  +  2  Q2 


(7) 


and  whether  or  not  ym  exists,  we  can  use  the  right  hand  side  of  (7)  (which  is  always  defined)  to  introduce 

,  ...  Vm 

k  =  1/sin  (— )  . 


Notice  that  the  representation  for  C  in  terms  of  the  initial  conditions  implies  that 
MV  +  2<2 2( 1  -  c°sy0) 


s'n‘(— )  = 


4  Q2 


and  so  it  will  follow  that 
signU)  =  sign(y0)  . 
Hence, 


V  =  4Q“[sim  (— )  -  sin-  (^) 


4 Q~t,  .  t  y 
=  — [!  -  ^sin-(|) 
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.  an  separate  the  v  ariables  and  integrate 


.  y  '  '0 


It  'a c  separate  the  Integra!  ■  n  the  left  side  into  two  integrals  and  make  the  substitution  $ 

obtain  the  ti'llowme  relati'Mi: 


=  then  we 


2m  e 

k 


|*  ?  2d<|>  _  f 2  2d0 

'/l  -  A2sin20  0  Vl  -  /f2sin2<J 


Q 

u  =  -  r(1i 


,a.^, 

2  Vl  -  *2sin2 


w  here  F(k.  *)  denotes  an  incomplete  elliptic  integral  of  the  first  kind. 

The  solution  gets  a  little  complicated  at  this  point,  for  if  there  does  exist  a  maximum  deviation  from  the 
libration  focal  point,  then  |£[ >  1  and  we  must  use  the  transformation 
F(k.  0)  =  if  (1.  sin  _,(A:sin 


Hence,  in  this  situation  we  must  transform  (8)  to 

1  ,  Vo 

u  =  Q(t  -  /0)  +  F  sin  U'sin— )) 


(*siny) 


/ 1  -  —  sin20 
k1 


In  the  case  where  |A|  >  1,  we  will  say  that  y  is  in  libration.  and  when  \k\  <  1  (i.e.,  the  right  hand  side  of 
(7)  >  1).  then  we  will  say  that  is  in  circulation.  In  this  setting,  y  will  'go  over  the  top'.  If  \k\  -  l.then 
V  will  be  in  the  rare  situation  of  being  exactly  'stalled'  between  circulation  and  libration.  We  will  ignore 
this  possibility  since  it  is  highly  unlikely  that  it  will  happen  in  a  situation  involving  artificial  satellites 
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(and  even  if  such  a  situation  happens,  then  the  nongravitational  effects  will  quickly  force  the  object  into 
either  circulation  or  libration).  Garfinkel  [2,  p.  660]  has  pointed  out  that  in  this  case  the  elliptic  functions 
degenerate  into  hyperbolic  functions,  and  so  an  analytical  solution  is  possible  in  this  case  as  well. 

In  the  case  that  |A  |  <  1,  then  there  will  not  exist  and  V) f  will  be  allowed  to  circulate  from  0  to  2k. 
The  complete  solution  for  V  can  be  obtained  from  (8) 

sin  (^)  =  sn(/<,  k)  ,  (10) 

where  sn(w.  k)  is  a  Jacobi  elliptic  function.  There  is  a  large  amount  of  literature  on  the  elliptic  integrals 
and  elliptic  functions,  and  calculation  to  any  desired  accuracy  is  possible  by  quickly  converging  series.  A 
good  reference  to  consult  for  the  complete  theory  of  these  transcendental  functions  is  Whittaker  and 
Watson  [8],  In  the  libration  setting,  the  complete  solution  of  \|/  can  be  obtained  from  (9)  and  is 


7 


3.  THE  SOLUTION  OF  THE  LAGRANGE  EQUATIONS  OF  MOTION 


If  we  become  a  little  less  restrictive  than  we  were  in  the  introduction  and  allow  q  to  be  not  zero  with  0) 

allowed  to  vary  (but  having  no  acceleration),  and  if  we  isolate  ourselves  to  one  critical  tesseral  harmonic 

in  the  disturbing  function,  then  we  establish  the  following  differential  equation  similar  to  (4): 

vji  =  ULm  «  pf  sin  y  for  /  -  m  even 
Sq  L-cosy  for  /  -  m  odd 

where 

and  we  have  dropped  the  subscripts  of  vy.  It  is  the  product  F  .(i)G  tnq(e)  which  determines  the  sign  of  P, 
and  so  by  adding  a  multiple  of  n  or  of  n/2  to  V,  if  necessary,  we  can  set 


and  produce  the  form  y  =  -0-sin  V  as  in  (6).  Our  assumptions  insure  that  Q  is  a  positive  constant,  and 
so  the  solution  outlined  in  the  previous  section  is  applicable.  It  is  interesting  to  note  that  the  application 
of  the  simple  pendulum  to  resonance  problems  dates  back  to  Laplace,  who  in  1800  was  investigating  the 
Galilean  satellites. 

The  Lagrange  equations  of  the  general  perturbations  of  the  Keplerian  elements  are  listed  in  Kaula 

[4.  p.  29)  and  are  reproduced  here  for  the  convenience  of  the  reader. 

da  _  2_  9K 
dr  na  dM 

de  _  1  -  e2  dR  "'ll  -  e2  dR 
nea ^  dM  nea^  dco 


dco 

d7 


dR 

de 
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dr 

cos  i 

dR 

1  dR 

It 

na24 1  -  e2sin  i 

dco 

na‘ 

Ml  -  e2sin i  dQ 

dfl 

1 

dR 

d It 

na2^!  1  -  e2sinr  ^ 1 

m 

1  -  e2 
=  n  -  - 

nea2 

dR 

2 

dR 

dr 

de 

na 

da 

The  n  appearing  on  the  right  hand  side  denotes  the  mean  motion  (expressed  in  the  proper  units).  It  is 

clear  that  if  we  remain  consistent  with  our  first  order  development  and  continue  to  hold  a,  e,  and  i 

constant  on  the  right  hand  side,  then  we  can  separate  the  above  system  of  differential  equations  and 

integrate  directly.  Keeping  in  mind  that  R  in  (3)  is  restricted  to  one  critical  tesseral  harmonic,  it  is  clear 

that  in  each  of  the  Lagrange  equations  we  must  either  evaluate  f '  sin  xj/dr  or  else  [ '  cos  ydr  .  We  will 

J'0  J'0 

proceed  to  demonstrate  how  this  is  done. 

Case  1  (|A|  >  1)  Libration. 


Recall  that  in  this  case  fcsin^  =  sn  (rr.i)  .  In  this  regime,  we  can  recenter  the  longitude,  if  necessary,  so 
that  ¥  is  restricted  to  the  first  and  fourth  quadrants,  and  so  cos  (^)  =  7 1  -  sin  2  (^)  .  With  the 
representation  of  sin  ^  j ,  an  identity  of  elliptic  functions  yields 

cos  (^)  =  VI  -  —  sn2 
2  k2 

Here  we  have  relied  on  the  fact  that  dn  (m.  I)  is  always  positive  when  the  argument  u  and  the  modulus  1 

are  real.  For  the  elementary  properties  of  elliptic  functions  that  will  be  assumed,  the  reader  should 
consult  [8]. 


(«.!)- Vdn’o.!)  =  dn  (m.  1) 


We  can  now  evaluate 


[  sin  \|/dr  =  [  2sin  A  cos  (^)d/  =  f  -sn  (n.  !)  dn  («, -)  dr 

J'o  J'o  2  2  V  k  k 

=*y,;n<"  j)dn<"' j,d"=©[,:n("»-r) -cnK)]  ■ 
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where  u  =  u(t)  and  w0  =  w(/0).  Furthermore, 

J*  cos\|/d/  =  J  2cos2|^j  -  Id/ =  J  2dn2^/,ij  -  Id/ 

- h^2 H)-ld-5l£K)-£ ("°- /)]  ~h (“ ' “»)  ■ 

where  £(«,  -!-)  is  the  fundamental  elliptic  integral  f“dn2ndw  [8,  pp.  517-518]. 
k  J  o 

Case  2  (|£|  <  1)  Circulation. 

Recall  that  in  this  case  sin  (^)  =  sr.(n,  it).  We  must  consider  u  in  (8)  as  a  uniformizing  variable.  As 

time  increases,  either  u  always  increases  or  always  decreases  (depending  upon  the  sign  of  k).  In  a  like 
manner,  V|/  will  always  increase  or  always  decrease  depending  upon  the  sign  of  \J/0  (which  is  the  sign  of  k). 
Hence  it  will  always  be  the  case  that  u  and  \j/  are  moving  in  the  same  direction  as  time  increases;  that  is, 
^  is  always  positive  in  the  circulation  regime.  Now  if  we  differentiate  both  sides  of  the  equation 

sin¥  =  sn(M.  k)  by  u.  we  observe  that 


cos 


and  since  dn(M,  k)  is  always  positive,  it  follows  that  sign  (cos^)  =  sign[cn(/<.  £)].  Moreover 
2 (^2)  =  *  ~  S*n2(^)  =  *  _  sn 2(u.k)  =  ctr(U'k)  , 


cos' 


and  so  it  follows  that  cos  (^ )  =  cn(«,  k). 


We  are  now  in  a  position  to  evaluate  the  integral  of  sinyd/: 


f  sin\|/d/  =  f  2sin  /' 

2 )  cos  ( ~ 

^]d/=  f  2sn(«.  k)cn{u,  k)dt 

J,0  J,0  ' ' 

L!  \  ■ 

l}  J'o 

2  k[u  2 

=  —I  sn(«,  k)cn(u.  k)du  =  — [dn(w0.  k)  -  dn(n,  k)] 


Qk 


li 


For  the  integral  of  cos  yd/,  we  begin  by  writing 

f  cos  yd/  =  f  1  -  2sin2(^d/  =  1  -  2sn2(«,  k)du  , 

J'0  J'o  '2'  QJ“0 

and  we  consider  the  identity  [8,  p.  516]  u  -  k^[usn2udu  -  fudn2ud«.  Hence, 

J0  Jo 

f  cos  yd/  =  !Uu  -  m0]  -  ^-|(w  -  u0)  ~  [£(“.  k)  -  E(u0 ,  k)}\ 

*  /q  Q 

=  jL[E(u,  k)  -  E(u0,  k)]  -  k^±liu  -  uQ)  , 

where  A'2  =  1  -  A-2  denotes  the  complementary  modulus. 

The  mean  motion,  n ,  on  the  right  hand  side  of  the  Lagrange  equation  for  the  rate  of  change  of  the  mean 
anomaly  should  not  be  assumed  constant  in  a  first  order  treatment;  for  the  mean  motion  appears  in  the 
equation  for  the  rate  of  change  of  M  in  the  unperturbed  situation  also.  Since  n  =  p1/2a_3/2,  where  p  is 
constant,  we  can  use  the  solution  for  a  outlined  above  to  compute  the  integral  V  ndt. 

In  our  restricted  resonance  problem,  we  have  that 


las, 


a  = 


Oy-,7 


3nm 


Q~ siny  , 


and  if  we  write 
4  as0Q 
3nmk 


A  = 


then  it  will  follow  that 


a  =  aQ  +  /\{cn  (“o-  j)  "  cn  (“•  for  1*1  >  1  . 

dn  ^uQ,  k'j  -  dn  («,  k j  for  |A|  <  1 

We  can  write  a  =  aQ  +  8fl,  and  use  the  binomial  series  to  approximate  n.  We  will  write  a  =  an  [l  +  —  1 

l  ao\ 

and  expand 


•  ■  -I1  - 1©  *  ¥©'  -  ]■ 


(12) 
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r 


from  which  point  we  will  truncate  the  rest  upon  the  assumption  that  the  series  drops  off  rapidly.  There 
are  then  four  integrals  that  we  need  to  evaluate,  and  with  reference  to  [8,  p.  516],  we  will  list  these 
integrals  evaluated: 

/„“("•  i)d“  =  [s"("i)'“"("'j)]  • 

ru 

I  dn(w,  k)du  =  am(«,  k)  =  sin-1  (sn(u.  A)}  +  j(u)n,  j(u)  integral  , 

*0 

where  k'z  =  1  -  — ,  and 

cu 

I  dn2(«,  k)du  =  E(u.  k)  -  £(«0,  k)  . 

J“0 

The  integral  f '  ndi  can  now  be  readily  evaluated  with  the  aid  of  the  above  integrals.  Hence,  each  of  the 
.  '0 

Lagrange  differential  equations  of  the  general  perturbations  of  the  Keplerian  elements  can  be  separated 
and  integrated  to  yield  the  complete  solution  of  the  restricted  resonance  problem  to  first  order. 

Suppose  that  the  disturbing  function  is  isolated  to  one  critical  tesseral  harmonic. 

R  =  ~~  (f)  ^Fimp^Glpq{e)[]lnfos  V  , 

where  y  =  m(k  -  \lm)  -  qi 0  +  jn  and  X  =  1  (M  +  co)  -  (6  -  Q).  Here  j  =  0.  1,  1/2.  or  3/2  depending 

*0 

upon  the  sign  of  F lmp(i)G tpq(e)  and  whether  I  -  m  is  even  or  odd.  Because  this  sign  is  important  in  the 
equations  which  follow,  we  will  use  a  symbol, 
a  =  sign {Flmp(i)Glpq(e)}  . 

Upon  calculating  the  rates  of  change  of  each  of  the  Keplerian  elements  given  by  the  Lagrange  equations 
of  motion  and  applying  the  integration  outlined  in  this  section,  we  arrive  at  the  explicit  expressions  for  the 
changes  of  each  of  the  Keplerian  elements  to  first  order.  For  the  convenience  of  the  reader,  we  will  list 
these  changes,  where  the  elements  n,  a,  e ,  and  /  appearing  on  the  right  hand  side  of  each  equation  are  held 
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fixed  and  constant  from  the  initial  osculating  elements.  To  save  writing,  we  will  allow  *  to  denote  the 
modulus  k  if  | Xr|  <  1  and  \/k  if  | k \  >  1  and  we  will  write  F',mp(i)  and  G' tpq(e)  to  denote  the  derivative  of 
these  functions  with  respect  to  their  argument. 

ba  =  a! cn(w0’  -  cn(“’  *)  for  1*1  >  1 

—  Hn (u  for  \Lr\  <  1 


Ldn(WQ.  *)  —  dn(«,  *)  for  |fc|  <  1  , 

2  f  cn(//0,  *)  -  cn(M,  *)  for  |&|  >  1 
*^^-dn(«Q,  *)  -  dn(u,  *)  for  |it|  <  1  , 


m  1  -  e2  -  Vl  -  e2  q^l  1  -  e2. 


^  ~  e2F,mn(i) 


fG'lpq(e) 


Vl  -  e2 


f  ^[£(m,  *)  -  E(u0,  *)]  -  1(m  -  m0)  for  |X:  |  >  1 
A[£(M<  *)  -  E(u0,  *)]  -  k  ^  \u  -  u0)  for  |X|  <  1  . 


&•«  Jl 

na3V 1  -  t 


( 

.  1  \1 

qcosi  +  m  1 

1  -  —cos  / 

X 

Jo 

i_f  cn(«0,  *)  -  cn(«,  *)  for  |X|  >  1 
*2  L(jn(M0,  *)  -  dn(«,  *)  for  |X  |  <  1 


^[£(m,  *)  -  £(m0,  *)]  -  1(m  -  m0)  for  |/t|  >  1 
-^[£(m,  *)  -  E(uq,  *)]  -  - ^  -■(»  -  m0)  for  |X|  <  1  , 


8M  =  J  (1  -  M  J cn(«0,  *)  +  Jcn2(M().  *)V,  _  t  ) 

(\  2"  *-dn («0,  *)  8a2  Ldn2(«0.  *)/ 
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+ 


&/3/4  _  15/4*-  Cn(«0,  *)\  F  tan  '[*sn(w.  *)/dn(w,  *)]  -  tan  '[*sn(M0.  *)/dn(w0,  *)] 

Q  \2a  4a2  L dn(«0*  *)/  *)  -  am(u0,  *) 

+  15fc42  f  *[£(« ,  *)  -  £(m0,  *)  -  E'\u  -  w0)] 

8 Qa2  *)  _  E(u0,  *) 

~  ^5(7)  ‘J,m  [«/  +  W,mp«)Glpq(e)\  -  o(l  -  e2)F/mp(/)(^^)]  x 

^[£(w,  *)  -  £(w0,  *)]  -  1(m  -  i/0)  for  |*|  >  1 
-^[£(u,  *)  -  £(u0<  *)]  -  - -  uQ)  for  |*|  <  1  . 

The  reader  can  observe  that  singularities  occur  in  the  above  equations  as  the  inclination  approaches  0  or 
as  the  eccentricity  approaches  0  within  the  precision  of  the  computer.  In  the  case  of  near  zero 
eccentricity,  since  the  eccentricity  function  G!p<j(e)  has  order  we  can  choose  q  =  0  so  that  the 
singularity  in  8e  is  removable.  Moreover,  when  q  -  0,  it  can  be  seen  from  Kaula's  expression  for  the 
eccentricity  function  [4,  p.  37]  that  G'/p0(e)  has  order  e,  and  so  the  case  of  the  singularity  with  zero 
eccentricity  is  of  no  consequence.  Finally,  in  the  case  of  near  zero  inclination,  the  singularity  can  often  be 
removed  with  the  aid  of  the  inclination  function  £/  (/  V  This  is  the  case  with  the  more  dominant  tesseral 
harmonics  such  asJ22  and-/31,  and  so  this  singularity  can  be  handled  by  avoiding  the  critical  harmonics 
for  which  the  inclination  function  does  not  remove  the  singularity  with  sin  / ,  /  =  0. 

It  was  noted  earlier  that  the  elliptic  functions  can  be  computed  by  rapidly  converging  series.  A  routine 
was  developed  for  computing  the  elliptic  functions  which  was  based  upon  the  method  outlined  in  the 
appendix  of  [6],  This  appendix  evidently  was  based  upon  a  section  of  Whittaker  and  Watson’s 
development  of  theta  functions  [8,  p.  486],  and  while  some  argue  that  it  is  convenient  to  expand  the 
functions  in  terms  of  the  modulus  *,  this  algorithm  transforms  the  modulus  to  a  smaller  parameter  q ,  and 
so  the  resulting  expansions  will  always  converge  rapidly  even  for  the  troublesome  case  where  the 
modulus  approaches  one.  This  method  was  able  to  compute  each  of  the  elliptic  functions  to  an  accuracy 
of  16  digits,  and  the  algorithm  was  extremely  fast.  The  function  am(M,  *),  which  up  until  this  point  has 
been  used  but  was  defined  in  a  somewhat  ambiguous  fashion,  will  be  expressed  more  precisely  in  the  next 
section  on  mean  elements. 
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4.  THE  MEAN  ELEMENTS 


In  this  section,  the  mean  elements  will  be  isolated  from  the  solution  outlined  in  the  last  section.  The 
periods  of  the  corrections  to  the  secular  changes  will  be  given,  and  it  will  shown  that  a,  e ,  and  i  all  have 
constant  secular  change,  while  the  angulai  arguments  to,  £2,  and  M  all  have  a  nonzero  secular  rale  of 
change.  It  should  be  noted  here  that  Garfinkel  [2]  expressed  his  solution  in  terms  of  mean  elements, 
which  he  obtained  by  integrating  the  elliptic  functions  over  their  periods  (as  we  will  do  shortly).  Also, 
Gedeon  [3]  gave  an  interesting  comparative  study  of  the  long  periods  among  satellites  with  varying 
parameters.  The  mathematics  involved  at  this  point  is  not  difficult  due  to  the  friendly  nature  of  the 
elliptic  functions  [8]. 

Let  us  introduce  the  number  K  =  F(*,  ml),  where  (as  before)  the  modulus  *  is  assumed  to  be  k  if 
|£|  <  1  and  1  Ik  if  |£|  >  1.  The  period  of  the  functions  sn(u,  *)  and  cn(u,  *)  is  4K,  whereas  an(u,  has  the 
smaller  period  2K.  Both  sn(w,  *)  and  cn(u,  *)  have  an  average  value  of  0  over  a  4K  period,  and  Garfinkel 
[2,  p.  663]  has  noted  that  the  average  value  of  dn(u.  *)  over  a  2K  period  is  n/2K.  Since  this  is  an 
important  fact  and  requires  the  evaluation  of  am(w,  *)  at  any  value  of  u,  we  will  present  a  proof  of  this 
fact. 


Recall  that  in  the  last  section  we  noted  that  f  Udnudu  =  am«.  Since  -nil  <  sin”1(snw)  <  nil,  we  must 

J0 

take  into  account  the  periodic  nature  of  sn  in  order  to  write  an  exact  formula  for  am(n,  *).  Let  r  denote 
the  integer  part  of  u/K  and  let  s  =  r(mod  4).  Then 


am(u,*)  = 


sin_I[sn(u,  *)]  +  ^(r  -  s)  if  s  =  0 

it  -  sin_l[sn(u,  *)]  +  H(r  -  s)  if  s  -  1  or  r  =  2 

2n  +  sin_1[sn(w,  *)]  +  5(r  -  s)  if  s  =  3  . 


Therefore  am(w,  *)  =  —  +  { am(u,  *)  -  — } .  where  the  part  in  braces  is  the  periodic  part  with  period  4K 
and  average  0,  and  nu/2K  is  the  secular  portion.  Hence, 


1  f2*  1 

—  dn(u,  *)d u  =  — am(2/t\  *) 


^  "  sin_l[sn(2/C,  *)]  -  ji] 
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We  also  require  the  secular  parts  of  £(«,  *)  and  tan  _1  [*sn(w,  *)/dn(«,  *)].  Now  E(u,  *)  is  not  a  periodic 
function,  but  it  can  be  expressed  by  the  relation  £(«,  *)  =  Z(m,  *)  +  Eu/K ,  where  Z(u ,  *)  denotes  the 
so-called  Zeta  function  (which  has  a  period  of  2K  and  average  0)  and  £  is  a  complete  elliptic  integral  of 
the  second  kind;  i.e.,  £  =  E(K,  *).  The  function  tan_1[*sn(w,  *)/dn(«,  *)]  is  periodic  with  periods,  and 
the  average  over  its  period  is  0,  although  this  fact  is  not  easy  to  prove. 

In  order  to  write  the  mean  elements  analytically,  we  must  have  analytical  expressions  for  <  [ f  sin  yd/>, 

J'o 

<f‘  cosyd/>,  and  <1 '  ndt>,  where  <•>  represents  the  part  of  •  with  the  periodic  portions  removed  by 
'0  J'o 

averaging.  In  order  to  write  these  explicitly,  we  must  again  divide  the  work  ir,.o  two  cases. 

Case  1  (|*|  >  1)  Libration. 


Since  v.Cii  (u,  -  )> 


-  G  over  a  4K  period,  we  have  that 


The  function  £(«)  has  a  nonzero  secular  rate  of  change,  and  so  for  the  average  of  the  integral  of  cosy,  we 
isolate  the  linear  pan  and  integrate  the  periodic  part  over  its  period  of  2K.  This  becomes 


The  secular  part  of  the  integral  of  the  mean  motion  n  is  more  lengthy  because  of  the  binomial  series 
expansion: 


18 


Case  2  (J  A'  |  >  1)  Circulation. 


Since  the  average  of  dn(u,  k)  over  a  2K  period  is  n/2K,  it  immediately  follows  that 


Similar  to  case  1,  we  have 


r!  7  i  \~>e 

"X cos V“'>  -  Tq  [it  £("»'  *>]  +  i \y  - (i 1  +  ‘ _> 


« -  w 


and  finally,  we  have 


\f  "d,>  =  ~  TA2Mu°-  A)]  “  am("°’  k) 

+  nQ  <|l  -  ^dn(«0,  k)  +  ip42dn2(t<0.  k) 

*  *> 


JL  +  -  ,„> 

2/f  8K  j  0 


We  are  now  able  to  observe  from  the  Lagrange  equations  for  the  general  perturbations  of  the  Keplerian 
elements  that  in  each  of  the  two  cases,  the  secular  rate  of  change  is  0  for  a,  e,  and  /  and  nonzero  for  CO,  Q. 
and  M. 
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5.  THE  NUMERICAL  STUDY 

The  method  described  in  the  preceding  sections  was  used  in  a  FORTRAN  computer  program  to 
compare  the  long  period  propagation  of  the  Keplerian  elements  with  a  numerical  integration  of  the 
Lagrange  equations  of  motion  (with  which  the  resonance  effects  are  accounted  for  automatically).  A 
variable  step-size  Adams  predictor-corrector  polynomial  integrator  was  used  for  this  comparison,  and  the 
test  was  performed  on  a  Harris  800  machine  with  quadruple  precision  (96  bit  arithmetic).  The  integrator 
was  originally  coded  by  Fred  T.  Krough  at  the  Jet  Propulsion  Laboratory  in  Pasadena.  California  in  April 
of  1969,  and  it  was  adapted  for  Harris  computers  by  E.  Mike  Gaposchkin  at  MIT  Lincoln  Laboratory  in 
November  of  1982.  The  objectives  of  this  study  were: 

1 .  to  estimate  the  accuracy  of  the  analytical  theory  over  the  long  resonance  period, 

2.  to  investigate  the  parameters  which  effect  the  accuracy. 

3.  to  observe  the  order  of  magnitude  of  the  perturbation  due  to  a  critical  tesseral  harmonic, 

4.  to  determine  the  decay  in  accuracy  when  two  or  more  critical  tesseral  harmonics  are 
included  in  the  disturbing  function  (and  it  is  assumed  in  the  analytical  solution  that  no  two 
tesseral  harmonics  interact),  and 

5.  to  numerically  compare  the  perturbation  produced  by  one  dominant  critical  harmonic  with 
that  produced  by  more  than  one  critical  harmonic. 

5.1.  AN  ERROR  ANALYSIS 

Before  the  comparative  results  are  presented,  it  is  necessary  to  optimize  the  analytical  model  in  the 

sense  of  minimizing  the  absolute  error  of  the  analytical  solution  of  the  total  perturbation  from  the  same 

perturbation  computed  by  numerically  integrating  the  Lagrange  equations  of  motion.  In  the  analysis  of 

the  solution,  it  is  clear  that  one  has  complete  flexibility  in  the  choice  of  the  elements  a,  e,  and  i  which  are 

held  fixed  and  constant  in  the  solution.  Apart  from  the  use  of  \j/0  and  y0  in  the  statement  of  the  initial 

value  simple  pendulum  problem,  one  has  freedom  to  change  the  choice  of  a,  e,  and  i  in  the  calculation  of 

the  modulus  k,  the  factor  Q  in  the  definition  of  u  and  in  the  evaluation  of  the  integrals  f '  sinvydf  and 

J,0 
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f '  cos  yd/,  and  in  the  factors  of  sin  y  and  cosy  in  the  Lagrange  equations  for  the  rates  of  change  of  the 

% 

elements.  For  notation,  we  will  refer  to  a  generic  Keplerian  element  (except  the  mean  anomaly)  by  a  and 
write  e  to  denote  the  factor  of  sin  y  or  of  cos  y  in  the  Lagrange  equation  for  the  rate  of  change  of  a. 
Therefore,  we  can  write  a  =  esiny  or  a  =  ecosy,  and  one  holds  e  =  e(a,  e,  i)  constant  over  an  entire 
resonance  period  and  proceeds  to  integrate  sin  y  or  cos  y  by  computing  k,  Q,  and  the  elliptic  functions  as 
outlined  in  Section  3.  It  should  not  be  surprising  for  one  to  observe  that  the  mean  elements  can  improve 
the  model.  Indeed,  the  mean  elements  (a),  ( e ),  and  (/)  are  constant  over  an  entire  resonance  period,  and 
we  will  see  that  if  they  are  positioned  strategically  in  the  analytical  solution,  then  a  significant 
improvement  is  realized  in  a  comparison  with  the  numerical  solution. 

We  can  iterate  on  the  integration  procedure  by  computing  the  mean  values  (6 a(a,  e.  /)>.  (5 e(a,  e,  i)>.  and 
(8 i(a,  e ,  /))  as  outlined  in  Section  4  and  on  each  successive  iteration  compute  8 a((a),  ( e ).  (/)),  where 
( a)=a0  +  (8a(a ,  e.  /)),  (e)  =  e0  +  (8 e(a.  e.  /')),  and  (/')  =  /0  +  <8/(a,  e,  /))  and  the  elements  a,  e.  and  i  are  saved 
from  the  last  iteration.  Empirically,  this  procedure  converges  rapidly,  and  in  fact,  only  two  iterations  are 
essentially  needed  to  realize  the  benefit. 

Consider  the  error  in  the  approximation  of  8a  for  a  synchronous  object  in  libration  with  the  harmonic 
^2 200'  In  F'gure  5-1  we  have  plotted  the  error  first  for  the  situation  where  the  initial  osculating  elements 
a0,  eQ.  and  iQ  are  used  to  determine  k,  Q,  and  e  and  then  for  the  situation  where  (a),  (e).  and  (t)  are  used  to 
compute  k,  Q,  and  e.  One  can  observe  that  in  the  first  case  the  error  is  grow  ing  without  bound  and  in  the 
latter  case  the  error  is  not  only  bounded,  but  the  error  is  bounded  by  a  small  amount.  Thus,  in  the 
libration  situation,  the  use  of  the  mean  elements  optimizes  the  model  of  the  motion  in  a  pleasing  way. 

The  above  result  makes  sense  physically  since  the  mean  values  of  a,  e,  and  i  correspond  to  the  stable 
equilibrium  longitude  of  the  oscillation,  about  which  the  motion  is  completely  symmetrical. 

Unfortunately,  the  situation  is  not  as  nice  in  the  circulation  regime  since  no  such  point  of  symmetry 
exists.  In  fact,  it  has  been  observed  that  any  values  of  a.  e,  and  i  in  the  analytical  solution  will  cause  the 
error  to  grow  without  bound  over  time,  but  it  has  also  been  observed  that  it  is  possible  to  minimize  the 
growth  by  using  the  mean  values  in  strategic  places.  The  optimum  is  not  achieved  by  simply  using  the 
mean  values  in  k,  Q.  and  e  as  in  the  libration  situation,  but  a  study  of  seven  synchronous  objects  in 
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ERROR  RESIDUAL  FOR  THE  SEMIMAJOR  AXIS 


AN  OBJECT  IN  LIBRATION  WITH  V2200 


Figure  5-1 :  The  error  cttrx  e  w  ith  initial  elements  inpur  and  then  with 
the  mean  elements  used  to  compute  k.  Q.  and  £ 

circulation  with  the  harmonic  has  revealed  that  there  are  two  cases  where  the  savings  can  be 

achieved.  It  was  observed  that  if  the  value  of  ( [ '  sin  \j/dr)  is  positive,  then  the  smallest  growth  is 

J,o 

achieved  by  using  the  initial  osculating  elements  for  the  computation  of  k  and  Q  and  the  mean  elements 

for  the  computation  of  £,  and  if  <  f '  sin  \\tdt)  is  negative,  then  the  smallest  growth  is  obtained  by  using 

% 

the  initial  elements  for  the  computation  of  Q  and  the  mean  elements  for  the  computation  of  k  and  £.  This 
phenomenon  is  not  well  understood,  but  it  appears  to  be  consistent  w  ith  every  circulation  object  studied. 
By  using  the  change  in  semimajor  axis  6 a  again  as  an  example,  we  demonstrate  the  first  case  in 
Figure  5-2  and  the  second  case  in  Figure  5-3.  Incidentally,  the  object  in  Figure  5-2  is  very  close  to  the 
border  between  circulation  and  libration,  and  therefore,  the  error  curve  resulting  from  the  mean  input 
appears  bounded  only  because  the  original  is  bounded. 


5.2.  THE  ALGORITHM 

Armed  with  the  above  error  analysis,  we  will  outline  the  algorithm  for  computing  the  changes  in  each 
of  the  Keplerian  elements  due  to  an  isolated  resonance  harmonic.  This  algorithm  is  designed  to  minimize 
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ERROR  RESIDUAL  FOR  THE  SEMIMAJOR  AXIS 


AN  08JECT  IN  CIRCULATION  WITH  V2200 


Figure  5-2:  The  error  curve  with  the  initial  elements  input  and  then  with 
the  mean  elements  used  to  compute  £ 


ERROR  RESIOUAL  FOR  THE  SEMIMAJOR  AXIS 


AN  OBJECT  IN  CIRCULATION  WITH  V2200 


Figure  5-3:  The  error  cur\  e  with  the  initial  elements  input  and  then  with 
the  mean  elements  used  to  compute  k  and  £ 
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the  variance  from  the  numerical  integration  of  the  equations  of  motion,  and  tnis  algorithm  was  used  in  the 
numerical  study  which  follows. 

1 .  The  osculating  elements  a0,  e0.  i(J,  vy0,  and  are  assumed  given  with  an  epoch  time  tQ. 

2.  The  dependent  variables  k,  Q,  and  e  are  computed  as  functions  of  the  initial  inputs. 

3.  The  mean  elements  (a).  { e ).  and  (/)  are  computed  as  defined  above. 

4.  If  |  A.  |  >  1 ,  then  k.  Q.  and  £  are  to  be  recalculated  by  using  (a),  ( e ),  (/),  \|/0.  and  Vo- 

5.  If  |  k\  <  1.  then  if  <fr  sinvydr)  <  0,  then  k  and  e  are  to  be  recalculated  using  (a),  ( e ).  (/'),  \| /Q, 

J'0 


6.  For  the  change  in  the  mean  anomaly  5Af,  the  integral  [ '  ndr  is  to  be  computed  as  a  function 

J'0 

of  k.  Q.  ,4  (as  a  function  of  mean  elements),  and  the  initial  semimajor  axis  aQ.  Recall  that 
the  elements  a  and  n  appearing  are  initial  elements  because  of  the  binomial  series  expansion 
which  was  used. 

7.  The  variable  u  and  the  elliptic  functions  are  to  be  computed  so  that  each  of  the  changes  can 
be  computed  as  listed  at  the  end  of  Section  3. 

5.3.  THE  RESULTS 

To  realize  the  objectives  of  our  study,  four  satellites  were  selected,  three  of  the  satellites  were  in  a  near 
synchronous  orbit,  and  the  other  was  in  a  near  half-synchronous  orbit.  Each  of  the  satellites  selected  was 
in  a  deep  resonance  with  at  least  one  critical  tesseral  harmonic;  for  it  is  in  this  regime  w  h ’re  the  error  is 
more  inclined  to  be  significant.  The  objects'  parameters  were  taken  from  an  object  file  maintained  by 
Group  91  at  Lincoln  Laboratory.  The  objects  are  catalogued  by  the  U.S.  Space  Command  catalogue 
system,  and  therefore,  each  object  will  be  referred  to  by  its  catalogue  number  in  that  file.  The  objects 
along  with  their  (osculating)  orbit  parameters  at  the  epoch  cited  are  listed  in  Table  5-1. 

For  each  object  three  plots  will  be  presented.  It  was  observed  in  each  of  the  situations  we  considered 
that,  although  each  ot  the  Keplerian  elements  was  affected  by  a  critical  tesseral  harmonic,  only  the  mean 
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TABLE  5-1 


Object  Parameters 


Object  No. 

14867 

15181 

13636 

16885 

a  km 

42170.5898 

42161.7406 

42166.032 

26553.963 

n  rev/day 

1.0025 

1 .0028 

1 .00267 

2.00635 

e 

2.71  X  10-3 

1.961  X  10-3 

5.714  X  10-4 

0.741 

<u  deg 

348.875 

180.467 

350.703 

288.15 

i  deg 

1.597 

1.087 

1.816 

63.257 

n  deg 

85.081 

84.648 

104.407 

79.338 

M  deg 

236.463 

179.122 

306.277 

25.459 

k  deg 

73.778 

116.064 

345.24 

35.459 

epoch  yr,  day 

'87,  140 

'87,  138.259 

'87.  139.5 

'87,  139.9 

anomaly  M  and  the  semimajor  axis  a  exhibited  large  changes.  We  have  therefore  restricted  our 
discussion  to  the  perturbations  of  these  two  elements.  For  each  object  we  have  plotted  the  following: 

1.  the  absolute  error,  or  residual  (the  perturbation  obtained  analytically  subtracted  from  the 
perturbation  obtained  numerically),  for  the  semimajor  axis  and  the  mean  anomaly,  of  the 
restricted  resonance  problem  where  the  disturbing  function  is  isolated  to  one  critical  tesseral 
harmonic,  and 

2.  the  comparison  of  the  perturbation  in  the  semimajor  axis  for  the  situation  when  the 
disturbing  function  is  isolated  to  one  dominant  critical  tesseral  harmonic  (V^OO  f°r 
example)  with  that  for  the  situation  where  the  disturbing  function  includes  more  than  one 
critical  tesseral  harmonic. 

The  last  type  of  plot  is  generated  entirely  from  numerical  data  and  is  intended  to  show  not  only  the 
magnitude  and  amplitude  of  the  perturbations,  but  also  to  what  extent  a  dominant  critical  tesseral 
harmonic  (if  such  a  harmonic  exists)  is  able  to  approximate  the  global  problem. 


26 


Tine  three  near  synchronous  objects  were  selected  on  the  basis  of  the  perturbation  produced  by  the 
critical  tesseral  harmonic  V'jjoo-  Object  No.  14867  is  closely  commensurate  with  the  rotation  rate  of  the 
earth  and  is  in  a  secure  libration  about  a  stable  equilibrium  point.  Object  No.  15181  is  also  in  libration 
about  a  stable  equilibrium  point,  but  because  of  its  initial  longitude  position,  it  is  in  a  long  libration  about 
that  point  with  a  large  amplitude.  Object  No.  13636  is  presently  in  circulation,  but  it  is  so  close  to  the 
border  between  circulation  and  libration  that  the  presence  of  other  geopotential  harmonics  or  other  effects 
in  the  disturbing  function  (such  as  lunar  and  solar  effects  or  solar  radiation  pressure)  could  cause  it  to 
librate  about  the  equilibrium  point. 

Among  the  synchronous  objects  (No.  15181,  No.  14867,  and  No.  13636),  the  harmonic  V^qo  was 
considered  to  be  dominant,  and  the  other  harmonics  included  in  the  study  were  V3J  J0,  V,,™,  1/42)q,  and 
^4400-  Notice  that  q  is  zero  in  each  of  the  harmonics  listed  above  because  all  of  the  synchronous  satellites 
in  our  study  have  small  eccentricity.  It  should  also  be  noted  that,  with  each  harmonic,  the  singularities 
involving  small  eccentricities  and  small  inclinations  were  removable.  For  the  Molniya  satellite 
(No.  16885)  with  high  eccentricity,  the  harmonic  \^220— i  was  isolated  as  dominant,  but  it  will  be  shown 
that  the  harmonics  V2211  anc*  ^3210  a*s0  Pr°duce  effects  of  large  amplitude.  The  harmonics  I^q  and 
^  421-1  were  also  usec*  *n  the  °f tflis  object. 

The  local  (or  isolated)  residuals  will  be  presented  first.  In  Figures  5-4  to  5- 1 1 .  the  residuals  for  the 
same  object,  element,  and  epoch  are  combined  on  the  same  plot.  One  can  see  from  those  plots  that,  for 
the  most  part,  the  analytical  theory  has  performed  satisfactorily  in  the  isolated  resonance  problem  (for 
which  it  was  designed),  with  the  exception  of  the  residuals  for  No.  13636  due  to  the  harmonic  V2200  and 
for  No.  16885  due  to  the  harmonic  ^32 1  o-  of  these  objects  are  in  circulation  with  the  respective 

.armonic.  In  order  to  understand  the  nature  of  the  residuals  a  little  better,  we  will  consider  additional 
parameters  which  are  indigenous  to  the  analytical  solution.  For  each  object,  we  will  list  in  Table  5-2  the 
values  k,  1  Ik,  and  X  for  the  force  which  produced  the  largest  amplitude  in  the  residuals.  For  the  three 
synchronous  objects,  this  harmonic  is  V'2200,  and  for  Object  No.  16885,  this  is  the  harmonic 

Consider  the  residuals  for  the  synchronous  objects  (Figures  5-4  to  5-9).  Object  No.  14867  is  librating 
about  the  stable  equilibrium  point  75°E,  and  its  longitude  position  and  velocity  indicate  that  it  has 
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RESIDUAL  OF  THE  SEMIMAJOR  AXIS 


OBJECT  No  14867 


(Thousands)  TIME  (d) 

Figure  5-4:  The  isolated  resonance  problem  for  a  synchronous  object 
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RESIDUAL  OF  THE  SEMIMAJOR  AXIS 

OBJECT  No.  16181 
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Figure  5-6:  The  isolated  resonance  problem  for  a  synchronous  object 


RESIDUAL  OF  THE  MEAN  ANOMALY 
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Figure  5-7:  The  isolated  resonance  problem  for  a  synchronous  object 
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RESIDUAL  OF  THE  SEMIMAJOR  AXIS 


OBJECT  No  13636 


Figure  5-8:  The  isolated  resonance  problem  for  a  synchronous  object 


RESIDUAL  OF  THE  MEAN  ANOMALY 


OBJECT  No  13636 


Figure  5-9:  The  isolated  resonance  problem  for  a  synchronous  object 
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RESIDUAL  OF  THE  SEMIMAJOR  AXIS 

OBJECT  No.  16885 


(Thousands)  TIME  (d) 

Figure  5-10:  The  isolated  resonance  problem  for  a  semisynchronous  object 


RESIDUAL  OF  THE  MEAN  ANOMALY 


(Thousands)  TIME  (d) 

Figure  5-11:  The  isolated  resonance  problem  for  a  semisynchronous  object 


31 


TABLE  5-2 

Parameters  Associated  with  the  Analytical  Solution 

Object  No. 

14867 

15181 

13636 

16885 

k 

-5.295 

1.51449 

-0.9986 

0.9263 

1/k 

-0.1889 

0.66 

-1.0014 

1 .0796 

\0  deg/day 

-8.267  X  10-2 

3.12  X  10-2 

-2.361  X  10-2 

0.159 

Harmonic 

2200 

2200 

2200 

3210 

recently  passed  this  point  and  is  now  decelerating  west.  The  larger  residuals  for  No.  15181  and 
No.  13636  from  the  harmonic  l^oo  correspond  to  the  larger  amplitudes  in  the  perturbations  which  result 
from  their  respective  longitude  positions  and  velocities.  Object  No.  15181  is  currently  about  41°  east  of 
the  stable  node,  and  although  it  is  drifting  eastward,  it  will  soon  be  drifting  west  towards  that  point. 
Because  of  the  distance  it  has  to  travel,  the  energy  it  will  attain  when  it  reaches  the  stable  point  will  cause 
the  amplitudes  of  the  perturbations  to  be  considerably  more  than  those  of  Object  No.  14867.  Since  Object 
No.  13636  is  drifting  west,  is  in  circulation,  and  is  near  the  unstable  equilibrium  point  345°E,  its  energy 
will  also  accelerate  so  that  again  large  amplitudes  in  the  perturbations  will  result.  This  analysis  places  the 
residuals  for  the  three  synchronous  objects  in  an  appropriate  context  (the  residuals  for  No.  14867  are 
smaller,  but  so  is  the  perturbation).  For  the  semimajor  axis,  the  difference  in  the  amplitude  of  the 
perturbation  for  these  objects  can  be  seen  in  Figures  5-12  to  5-14. 

The  amplitudes  of  the  residuals  displayed  in  Figures  5-10  and  5-1 1  for  Object  No.  16885  due  to  the 

harmonic  are  larger  than  expected,  and  it  is  not  clear  why  these  are  much  larger  than  the  residuals 

due  to  the  harmonics  V220_i  and  V22|(.  The  object  is  in  circulation  with  the  harmonic  Vj^q,  has  a  large 

drift  rate,  and  the  resonance  period  is  much  longer  in  this  case  than  those  periods  from  V'220_,  and  V22 1 1 - 

It  is  possible  to  blame  this  on  modeling  error  since  the  object  has  high  eccentricity  and  inclination,  and  so 

•  •  •  • 

the  fundamental  assumption  A.  =  M/s0  may  not  be  valid  over  this  long  resonance  period.  On  a  more 
positive  note,  it  is  obsers  ed  that  the  effects  due  to  the  harmonics  causing  a  libration  (namely  ^220-1  and 
V22n)  appear  to  be  modeled  very  well  by  the  first  order  theory. 
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OBJECT  No.  14867  CHANGE  IN  SEMIMAJOR  AXIS 


V2200,  V3110.  V3300.  V4210.  AND  V4400 


Figure  5-12:  Perturbation  comparison  of  the  local  and  global  problem 


OBJECT  No.  15181  CHANGE  IN  SEMIMAJOR  AXIS 


(Thousands)  TIME  (d) 

Figure  5-13:  Perturbation  comparison  of  the  local  and  global  problem 
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OBJECT  No.  13636  CHANGE  IN  SEMIMAJOR  AXIS 

V2200.  V3110.  V3300,  V4210,  AND  V4400 


Figure  5-14:  Perturbation  comparison  of  the  local  and  global  problem 

Consider  now  the  plots  in  Figures  5-12  to  5-16.  These  graphs  are  generated  entirely  from  data 
computed  from  a  numerical  integration  of  the  equations  of  motion,  and  each  Figure  depicts  the 
perturbation  in  the  semimajor  axis  in  the  two  cases  of  the  disturbing  function  containing  one  critical 
tesseral  harmonic  and  of  the  disturbing  function  containing  five  critical  tesseral  harmonics.  Hence,  the 
interaction  of  the  resonance  effects  is  displayed. 

Of  the  three  synchronous  objects,  one  can  see  that  the  perturbation  from  the  global  resonance  problem 
is  best  approximated  by  the  isolated  harmonic  V2200  *n  t*ie  case  Object  No.  14867.  The  period  and 
phase  differ  to  a  small  degree,  but  the  amplitude  seems  to  match  up  exactly.  This  is  also  somewhat  true  in 
the  case  of  Object  No.  15181,  but  the  amplitude  is  several  kilometers  off  in  the  isolated  case  and  the 
period  is  different.  Still,  one  can  make  a  case  with  each  of  these  two  objects  that  the  isolated  resonance 
problem  describes  the  perturbation  reasonably  well  for  a  couple  of  hundred  days.  This  is  far  from  true  in 
the  case  of  Object  No.  13636.  Indeed,  the  global  perturbation  looks  almost  as  if  a  long  and  slow  libration 
about  the  stable  equilibrium  point  is  modeled,  whereas  we  have  already  seen  that  this  object  is  in 
circulation  when  the  disturbing  function  is  isolated  to  one  tesseral  harmonic.  Of  course,  this  problem 
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OBJECT  No.  16885  CHANGE  IN  SEMIMAJOR  AXIS 

V220-1 .  V2211.  V3210,  V4410.  AND  V421-1 


(Thousands)  TIME  (d) 


Figure  5-15:  Perturbation  comparison  of  the  local  and  global  problem 


OBJECT  No  16885  CHANGE  IN  SEMIMAJOR  AXIS 

V220-1.  V2211.  V3210.  V4410,  AND  V421-1 


(Thousands)  TIME  (d) 


Figure  5-16:  Perturbation  comparison  of  the  local  and  global  problem 
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arises  from  the  fact  that  Object  No.  13636  is  close  to  the  border  between  circulation  and  libration  and  the 
other  effects  may  cause  the  object  to  be  perturbed  in  an  entirely  different  fashion  than  predicted  by  one 
harmonic.  Hence,  the  motion  close  to  the  border  between  circulation  and  libration  is  very  difficult  to 
predict,  even  though  the  object  is  in  close  commensurabiiity  with  the  earth’s  rotation  and  has  favorable 
parameters  such  as  small  eccentricity  and  inclination. 

Finally,  the  plots  in  Figures  5-15  and  5-16  show  the  complexity  involved  when  an  object  is  in 
resonance  and  has  high  eccentricity.  Both  of  the  forces  V/220—  i  and  ^22)1  have  s'm''ar  amplitudes,  but  the 
periods  are  vastly  different,  and  both  periods  are  far  from  the  period  exhibited  in  the  global  perturbation. 
Therefore,  one  has  little  hope  in  approximating  the  global  problem  by  using  one  isolated  critical  tesseral 
harmonic  for  this  situation. 

We  are  led  to  the  question  about  the  performance  of  the  analytical  approximation  to  the  global  problem 
by  assuming  that  the  critical  tesseral  harmonics  do  not  interact.  Unfortunately,  an  unsatisfactory  model  of 
the  problem  is  made  over  an  entire  resonance  period.  Often  the  residuals  exceed  the  actual  amplitude  of 
the  total  perturbation.  This  is  displayed  in  Figure  5-17,  where  the  residual  of  the  semimajor  axis  is 
presented.  The  only  positive  point  that  can  be  made  here  is  that  in  each  situation,  a  good  model  of  the 
problem  io  made  over  one  to  two  hundred  days,  and  sometimes  this  is  better  than  if  we  relied  on  only  one 
harmonic  to  dominate  and  model  the  global  problem  over  this  time  span.  We  conclude  that  the 
interaction  between  the  harmonics  is  significant  even  in  the  most  favorable  situations. 
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OBJECT  No  14867  SEMIMAJOR  AXIS 

V2200.  V3110.  V3300.  V4210,  AND  V4400 


(Thousands)  TIME  (d) 

Figure  5-17:  Residual  for  the  global  resonance  problem 
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6.  THE  RESONANCE  ADDITION  TO  THE  ANODE  MODEL 


The  ANalytical  Orbit  DEtermination  routine  (ANODE)  at  Millstone  Hill  is  an  entirely  analytical 
propagation  model  which  is  tied  to  a  least  squares  best  estimate  routine.  This  software  is  designed  to  fit 
data  in  near  real  time  for  the  maintenance  and  upkeep  of  the  satellite  object  element  sets  in  our  data  bases. 
ANODE  has  been  working  on  a  routine  basis  for  the  past  seven  to  eight  years,  but  often  the  quality  of  the 
element  sets  has  been  poor  for  synchronous  and  half-synchronous  satellites  because  of  poor  data  and  the 
lack  of  resonance  in  the  propagation  model.  The  original  decision  to  leave  resonance  out  of  the 
propagation  was  due  to  fears  that  resonance  would  slow  down  the  routine  and  also  because  storage  space 
was  a  fundamental  concern  at  that  time.  Beside  those  concerns  is  the  fact  that  at  that  time  the  overall 
quality  of  the  data  required  that  each  object  be  updated  every  three  to  seven  days.  Since  the  resonance 
correction  would  be  very  small  over  that  time  frame,  it  was  not  too  difficult  to  reaquire  objects  with  the 
instruments  that  were  used.  Nevertheless,  problems  such  as  mistagging  objects  or  losing  (low-priority) 
objects  arose  from  sparse  ‘angles-only '  optical  data,  poor  quality  data,  and/or  the  lack  of  quality  radar 
range  data.  While  better  data  w  ould  certainly  improve  this  aspect  of  the  problem,  there  is  also  a  need  to 
improve  the  ANODE  physical  model.  Moreover,  the  quality  of  data  now  available  is  much  better  than  ten 
years  ago.  and  so  there  can  be  a  longer  time  span  between  observation  for  much  of  the  catalogue.  This 
w  ill  be  important  as  an  increasing  number  of  objects  in  space  place  a  greater  strain  on  the  space 
surveillance  network,  and  so  it  is  necessary  to  improve  the  propagation  model  in  ANODE  to  better  system 
performance. 

The  resonance  corrections  described  in  Section  3  have  been  added  to  ANODE,  and  this  section  will 
outline  how  this  has  been  accomplished.  Before  giving  more  detail,  it  is  important  to  understand  how  the 
ANODE  model  works  without  the  corrections.  We  will  use  a  boldfaced  z  to  denote  a  vector  containing 
the  six  Keplerian  elements. 

The  propagation  in  ANODE  can  be  described  by  the  following  simple  stepwise  procedure.  First,  it  is 
assumed  that  mean  elements  z0  with  an  epoch  time  are  input,  and  the  secular  rates  z(z0)  are  calculated. 
The  partial  derivatives  of  these  rates  with  respect  to  the  input  mean  elements  z 0  are  also  calculated 


dz(z0) 
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where  the  formal  notation  that  is  displayed  has  been  adopted  from  [7],  Now,  for  propagation  to  a 
particular  time  r,  ANODE  propagates  the  mean  elements 

=  z0  +  i{!  -  /0)  +  -  t0)2X  [^]A  ’  (13) 

adds  the  long  periodic  corrections 

zL(l)  =  zm(t)  +  AL[zm(t)]  04) 

as  a  function  of  the  mean  elements,  and  adds  the  short  periodic  corrections  to  produce  the  osculating 
element  set 

*oscW  =  Z/.(,)  +  AS^L^  ■  (15) 

The  long  periodic  corrections  &L  arise  from  averaging  the  equations  of  motion  over  the  period  of  the 

moon  (28  days)  or  the  period  of  the  motion  of  the  argument  of  perigee,  whichever  is  shorter,  and  the  short 
periodic  corrections  come  from  averaging  the  equations  of  motion  over  the  period  of  the  satellite. 

Now,  in  most  deep  resonance  problems  (where  the  order  of  magnitude  of  the  perturbations  is 
significant),  the  period  of  the  resonance  motion  is  on  the  order  of  three  years  or  more.  This  is  much 
lorger  than  the  long  period  average  described  above,  and  therefore  it  is  decided  to  add  the  resonance 
corrections  to  the  secular  (or  mean)  elements.  The  resonance  corrections  are  included  by  the  follow  ing 
additional  steps.  Resonance  first  must  be  initialized  (as  outlined  in  the  first  five  steps  of  Section  5.2)  as  a 
function  of  the  osculating  element  set  z<w  (r0)  [as  computed  above  in  (15)]  and  a  critical  index  set 
representing  a  particular  tesseral  harmonic.  Then,  in  order  to  propagate  to  a  specific  time  I.  the  vector  of 
resonance  corrections  8z  (as  produced  at  the  end  of  Section  3)  is  added  to  the  mean  element  propagation 
(13) 

zm "(t)  =  zm(t)  +  8z  ,  (16) 

the  long  periodic  corrections  are  computed  as  in  ( 14) 

*z>>  =  *„;<'> +  ^[z„>] .  d7) 

and  the  short  periodic  corrections  are  added  to  produce  the  osculating  element  set  as  in  ( 15) 

z01r*('>  =  zZ>  +  %[**>]  •  08) 
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For  the  least  squares  fit  procedure  in  ANODE,  the  analytical  Jacobian  of  the  computed  pointing  C(t) 
with  respect  to  the  input  mean  elements  Zq  is  needed  [7],  The  function  C(t)  is  a  vector  containing  the 
topocentric  coordinates  (azimuth,  elevation,  range,  and  range  rate)  which  point  to  the  object  from  a 
particular  site.  Since  C(t)  is  merely  a  geometrical  transformation  from  the  inertial  position  computed 
from  the  osculating  element  set,  it  is  pointed  out  in  [7,  pp.9,  lOjthat  the  Jacobian  can  be  computed  by 
using  the  chain  rule 


dCU) 

'  dCU)  ' 

‘dzosA'y 

3z  «(')■ 

.  3zo 

.  5z0  . 

where  it  is  assumed  that  the  matrix 


(19) 


dzmU) 


is  the  identity  (since  it  is  a  function  of  periodic  corrections!  and  we  can  write 


aznl<f> 

dz0 


I 


+ 


<'-'0>  ■ 


where  Ib  is  the  6x6  identity  matrix. 

After  the  resonance  correction  is  added,  we  need  to  change  ( 19)  to 


(20) 


(21) 


(22) 
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The  first  matrix  on  the  right  hand  side  of  (22)  is  computed  as  in  (20)  and  the  second  matrix  is  adequate  to 
approximate 


because  the  error  in  so  doing  is,  again,  a  function  of  periodic  corrections. 

Most  of  the  work  involved  in  computing  the  second  matrix  on  the  right  hand  side  of  (22)  is  of  a  routine 
nature  using  the  chain  rule  and  product  rules  (and  we  will  leave  this  to  the  reader),  but  for  the 
convenience  of  the  reader,  the  partials  of  those  elliptic  functions  which  are  relevant  to  the  resonance 
corrections  are  listed  below  in  Table  6-1. 


TABLE  6-1 

Partial  Derivatives  of  the  Elliptic  Functions 

Function 

d 

<9u 

d 

<9k 

sn(u,  k) 

cn(u,  k)  dn(u,  kj 

-cn(u,  k)  dn(u,  k)  A 

cn(u,  k) 

-sn(u,  k)  dn(u,  k) 

sn(u,  k)  dn(u,  k)  A 

dn(u,  k) 

-k2sn(u,  k)  cn(u,  k) 

dn(u  k)  ^sn^u'  **)  cn(u,  k)  dn(u,  k)  A  -  sn2(u,  k)] 

E(u,  k) 

dn2(u,  k) 

~  [E,(k,  d»)  -  F(k,  </>)]  -  dn2(u,  k)  A 

am(u,  k) 

dn(u,  k) 

-dn(u,  k)  A 

A  few  preliminary  comments  are  needed  before  the  table  of  partial  derivatives  can  be  presented. 
Suppose  the  argument  u  and  the  modulus  k  are  given  and  set 
$  =  am(n,  k) 

with 


24 
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and 


rO 

not.  <t>)  =  - 

J0(l 


-  /t2sin 


Since  the  integrand  in  the  definition  of  F  is  continous,  a  direct  computation  (using  the  fundamental 
theorem  of  calculus)  yields 


ffl*-**  =  (1  -  A'2sin24>)~,/2 
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and 


A  =  =  iincA.  <!»  -  F(k,  0)]  . 

ok  k 

Finally,  let  the  symbol  Ej(k.  0)  denote  the  incomplete  elliptic  integral  of  the  second  kind 

f  0 - - 

E/(k,  0)  =  \1  -  k-sm-'^A't,  . 

J0 

The  above  procedure  can  be  used  for  as  many  critical  tesseral  harmonics  as  desired  by  simply 
calculating  each  set  of  corrections  separately  and  adding  them  all  together,  but  since  no  interaction  is 
modeled,  the  routine  will  decay  in  accuracy  over  time  as  discussed  in  Section  5.  Nevertheless,  for 
surveillance  purposes,  w here  theory'  is  coupled  with  real  data  and  differential  correction,  the  routine 
should  be  trustworthy  within  100  or  200  days  from  the  epoch.  Since  the  other  aspects  of  the  ANODE 
propagation  model  are  not  as  accurate  for  this  long,  this  should  not  be  too  discouraging. 

As  an  example  of  how  the  resonance  addition  has  improved  ANODE,  a  simple  test  has  been  performed. 
For  an  object  in  a  synchronous  orbit,  we  generated  20  days  worth  of  simulated  metric  data  using  the 
precision  numerical  integration  program  DYNAMO.  Since  the  data  created  is  essentially  free  of  noise, 
we  were  able  to  ignore  the  problem  of  bad  data  so  that  we  could  concentrate  on  the  propagation  model  in 
ANODE.  For  a  force  model  in  DYNAMO,  we  limited  the  effects  to  those  caused  by  the  geopotential 
expansion  through  degree  /  =  4  and  order  m  =  4,  and  this  includes  all  values  of  p  and  q  for  each  pair 
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(/,  m).  Hence,  we  only  modeled  the  geopotential  effects  on  the  object  which  ANODE  modeled,  and  we 
are  deferring  the  effects  caused  by  the  moon  and  the  sun  to  a  later  study.  The  data  was  sampled  every 
quarter  of  a  day  so  that  a  thorough  representation  of  the  orbit  over  the  20  day  period  is  gathered.  As  a 
test,  using  the  fit  routine,  we  fit  this  data  to  the  ANODE  propagation  model  for  three  different  cases:  I  —  a 
geopotential  model  consisting  of  ^20-^30.  and  ^40-  II  --  a  geopotential  model  consisting  of  '/20*  ^30’  V40’ 
and  one  resonance  term  V22oo> HI  "  a  geopotential  model  consisting  of  V20,  V^q,  V4q,  and  three 
resonance  terms  ^2200’  ^31 10’  ^3300-  Case  I  is  the  old  ANODE  model  (without  the  lunar  and  solar 

effects),  case  II  contains  the  dominant  resonance  effect,  and  case  III  attempts  to  model  the  complete 
problem  (as  much  as  is  possible  analytically). 

Case  I  did  not  fit  the  data  very  well.  After  10  least  squares  fit  iterations,  much  of  the  data  near  the 
extreme  times  was  edited  out  and  the  root  mean  square  (rms)  residual  for  the  azimuth  was  0.045°  for  all 
the  data  and  0.017°  for  the  edited  data,  for  the  elevation  was  0.01 8°  for  both  sets  of  data,  for  the  range 
was  1.068  km  for  all  the  data  and  0.3878  km  for  the  edited  data,  and  for  the  range  rate  was  0.0047  m/s  for 
both  sets  of  data.  All  of  the  data  fit  nicely  for  case  II  and  case  III.  and  in  each  case  only  four  iterations 
were  required.  The  rms  residuals  for  case  II  were  0.004°  for  the  azimuth.  0.002°  for  the  elevation, 

0.0907  km  for  the  range,  and  0.0034  m/s  for  the  range  rate,  and  the  rms  residuals  for  case  III  were  0.001  ° 
for  both  the  azimuth  and  the  elevation.  0.0049  km  for  the  range,  and  0.0033  m/s  for  the  range  rate. 

Hence,  the  latter  model  is  able  to  recover  the  effects  to  within  a  millidegree  in  azimuth  and  elevation  and 
within  5  meters  in  range. 


7.  SUMMARY 


A  first  order  analytical  theory  of  resonance  effects  on  satellite  orbits  has  been  developed  and  tested. 
The  motion  is  modeled  by  the  motion  of  a  simple  pendulum  with  initial  conditions,  and  is  thus  naturally 
defined  in  terms  of  the  elliptic  functions  of  Jacobi.  The  theory  has  been  tested  against  a  numerical 
integration  of  the  equations  of  motion  for  the  general  perturbation  of  the  Keplerian  elements,  and  the 
following  conclusions  can  be  drawn. 

1.  The  theory  is  easy  to  manipulate  directly  to  produce  the  periods,  the  mean  elements 
averaged  over  the  resonance  periods,  the  periodic  corrections,  and  the  partial  derivatives  of 
the  resonance  corrections  with  respect  to  the  input  elements. 

2.  The  method  works  best  when  one  isolates  the  resonance  problem  to  that  of  a  single  critical 
tesseral  harmonic,  although  this  is  seldom  adequate  to  represent  the  ‘real  world’  problem. 

In  fact,  there  is  not  guaranteed  to  be  a  dominant  critical  tesseral  harmonic  even  in  the  case 
of  a  circular  equatorial  orbit  which  is  very'  near  synchronous. 

3.  The  global  problem  can  be  modeled  analytically  by  solving  the  isolated  resonance  problem 
locally  and  adding  the  effects  together,  but  this  is  not  a  good  model  over  a  full  resonance 
period,  and  is  only  able  to  describe  the  motion  w  ith  reasonable  accuracy  over  a  short  time 
span. 

4.  This  method  is  good  for  a  first  approximation  to  the  solution  of  the  resonance  problem,  but 
the  theory  should  not  be  used  for  accurate  predictions  without  iteration  correction  by 
comparison  with  other  data  (such  as  observed  data  or  numerically  simulated  data).  Hence, 
some  applications  of  this  theory  are: 

a.  as  a  front  end  processor  for  a  numerical  integrator, 

b.  for  addition  to  an  entirely  analytical  theory  which  fits  observed  data  in  real  time  and 
attempts  to  correlate  targets  quickly  (i.e.,  this  algorithm  has  been  implemented  in  the 
analytical  orbit  determination  routine  ANODE  at  Lincoln  Laboratory), 


c.  for  gaining  insight  into  a  physical  problem  involving  resonance  such  as  a  station 
keeping  pioblem  or  a  problem  wmch  requires  a  general  notion  of  the  long  term 
effect  of  a  particular  Keplerian  element  from  a  resonance  force.  An  example  of  this 
application  is  the  initial  modeling  of  the  resonance  effects  which  were  important  for 
the  design  of  a  series  of  station  keeping  maneuvers  for  LES-8  [5], 

5.  An  object  which  is  commensurate  with  the  rotation  rate  of  the  earth  can  be  in  deep 
resonance  with  several  tesseral  harmonics,  and  the  interaction  often  makes  it  difficult  to 
analytically  model  the  motion  either  by  assuming  no  interaction  or  by  hoping  one  tesseral 
harmonic  will  suffice.  Two  significant  parameters  which  determine  the  resonance  motion 
are  the  initial  longitude  position  and  the  initial  longitude  rate  (along  with  the  direction  of  the 
longitude  motion).  The  problem  is  compounded  by  the  reality  that  an  object  may  be  in 
libration  with  one  critical  tesseral  harmonic  and  in  circulation  with  others.  Of  course,  either 
libration  or  circulation  will  dominate,  but  the  damage  to  the  hope  of  modeling  the  global 
problem  by  solving  the  isolated  problem  separately  and  adding  the  effects  together  is  clear. 
Moreover,  it  is  clear  from  the  magnitude  of  the  amplitudes  displayed  that  the  resonance 
effects  on  satellite  orbits  should  not  be  ignored  when  an  object  is  commensurate  with  the 
rotation  rate  of  the  earth. 
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